%%
% var_case1

%% setup
%clear all;close all;clc
setup_workspace;
try
s = RandStream('mt19937ar','Seed',999);
RandStream.setGlobalStream(s);
catch
disp('seed did not work')
end

%% create a folder
ps.fpath.cucd=pwd; % current folder
cd(ps.fpath.Case)
mkdir(ps.file_name) % make folder
cd(ps.fpath.cucd)

%% load data
% us data
us_row_new1; % load use data
%change the order of variables
tmp=db.data1;
db.data=tmp([1 2 3 6 7  4 8 5 ],:); % TFP first
tmp=db.data_name;
db.data_name=tmp([ 1 2 3 6 7  4 8 5]);


db.data=db.data*100;%*25; % normalization of data
%db.data=db.data(:,29:280); % cut nan (from 1955q1 to 2017q4)

ps.T=size(db.data,2); % data length
ps.N=size(db.data,1); % number of variables  

%--------------------------------------------------------------------------

%% set parameters

ps.p=2; % number of lags in VAR
ps.cons=1; % constant in VAR if 1 ?
%ps.OD=1; % order of MA representation computed
ps.prior=2.5; % Minnesota prior = 2
ps.nsave=1000; % draws saved
ps.nburn=100; % burn in
ps.ndraws=ps.nsave+ps.nburn; % 
ps.nprint=100; %
ps.stationary=1; % VAR stationary 0 (do nothing), 1 (eliminate non-stationary)


% variable to focus on in VAR (k-th variable)
%--------------------------------------------------------------------------
% ps.k=1; % k-th variable in VAR for MBC shock
%--------------------------------------------------------------------------
ps.nbin=1*1024; % number of bins for approximation of integral

ps.w1=2*pi/32; % w1: lower value
ps.w2=2*pi/6; % w2: upper value
ps.grid_min=0;%1e-6;
ps.grid_max=2*pi;
ps.ngrid=1*1024;% number of grids in [0,2*pi];
ps.grid=linspace(ps.grid_min,ps.grid_max,ps.ngrid); % equally spaced grids
ps.igrid=zeros(1,ps.ngrid); % 1 if the grid is [w1,w2], ow 0
for i=1:ps.ngrid
    if (ps.grid(i)>=ps.w1)&&(ps.grid(i)<=ps.w2)
        ps.igrid(i)=1;
    end
    if (ps.grid(i)>=2*pi-ps.w2)&&(ps.grid(i)<=2*pi-ps.w1)
        ps.igrid(i)=1; 
    end    
end

%--------------------------------------------------------------------------

%% run VAR
% Bayesian VAR
out=bvar_f_w(ps,db);

% standard VAR
%out=VAR_OLS(db.data,ps.p,ps.cons,ps.OD);

%--------------------------------------------------------------------------

%% compute MBC shock q for the k-th variable
% reorder ALPHA (N by N by p) N=# of variables, p=# of lags
out.ALPHA_draws2=nan(ps.N,ps.N,ps.p,ps.nsave);% preallocatoin
for i=1:ps.nsave
    for j=1:ps.p
    tmp=out.ALPHA_draws(2:end,:,i); % take out constants
    out.ALPHA_draws2(:,:,j,i)=transpose(tmp(1+(j-1)*ps.N:j*ps.N,:)); % transpose
    end
end

for ii=1:ps.N 
    
% variable to focus 
ps.k=ii;
ps.k_w=1; % one variable
    
% compute q for each draw
rs.q_draws{ii}=nan(ps.N,ps.nsave); % preallocation
rs.eig_vec_draws{ii}=nan(ps.N,ps.N,ps.nsave);
rs.eig_val_draws{ii}=nan(ps.N,ps.N,ps.nsave);
rs.D_all{ii}=nan(ps.ngrid,ps.nsave);
rs.qDq{ii}=nan(ps.ngrid,ps.nsave);

% set options
opt.w1=ps.w1; % frequency
opt.w2=ps.w2; % frequency
opt.k=ps.k; % variable to focus for MBC shock identification
opt.k_w=ps.k_w; % weight for each variable

%tic
for i=1:ps.nsave % draws
% input for compute_MBC_shock 
tmp_out.BB=out.ALPHA_draws2(:,:,:,i); % coefficients
tmp_out.SIGMA=out.SIGMA_draws(:,:,i); % covariance

% compute q
%tmp_rs=compute_MBC_shock(tmp_out,ps,opt); 
tmp_rs=compute_MBC_shock4(tmp_out,ps,opt);

% keep draws
rs.q_draws{ii}(:,i)=tmp_rs.q; 
rs.eig_vec_draws{ii}(:,:,i)=tmp_rs.eig_vec;
rs.eig_val_draws{ii}(:,:,i)=tmp_rs.eig_val;

% spectral density
rs.D_all{ii}(:,i)=tmp_rs.D_all_vec;
rs.qDq{ii}(:,i)=tmp_rs.qDq;

end
%toc
end

%--------------------------------------------------------------------------

%% save 

ps.fpath.cucd=pwd; % current folder
cd(ps.fpath.Case)
mkdir(ps.file_name) % make folder
cd([ps.fpath.Case '/' ps.file_name]);
%save([ps.file_name '_var'],'allvar','-v7.3')
save([ps.file_name '_var_result'],'ps','rs','out','db','opt','-v7.3')

cd(ps.fpath.cucd)










